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Abstract 

In this paper, we consider two stochastic models of gene expression in 
prokaryotic cells. In the first model, sixteen biochemical reactions involved 
in transcription, translation and transcriptional regulation in the presence 
of inducer molecules are considered. The time evolution of the number of 
biomolecules of a particular type is determined using the stochastic simula- 
tion method based on the Gillespie Algorithm. The results obtained show that 
if the number of inducer molecules, Nj, is greater than or equal to the number 
of regulatory molecules, Nr, the average protein level is high in the steady state 
(state 2). The magnitude of the level is the same as long as iVj > Nr. When 
Ni is -C Nr, the average protein level is low, practically zero (state 1). As Ni 
increases, the protein level continues to remain low. When Ni becomes close to 
Nr, protein levels in the steady state are intermediate between high and low. 
In the presence of autocatalysis, a cell mostly exists in either state 1 or state 2 
giving rise to a bimodal distribution in the protein levels in an ensemble of cells. 
This corresponds to the "all or none" phenomenon observed in experiments. In 
the second model, the inducer molecules are not considered explicitly. An ex- 
haustive simulation over the parameter space of the model shows that there 
are three major patterns of gene expression, Type A, Type B and Type C. 
The effect of varying the cellular parameters on the patterns, in particular, the 
transition from one type of pattern to another, is studied. Type A and Type 
B patterns have been observed in experiments. Simple mathematical models of 
transcriptional regulation predict Type C pattern of gene expression in certain 
parameter regimes. The model studied by us includes all the major biochemi- 
cal reactions involved in gene expression and the stochastic simulation results 
provide an understanding of the microscopic origins of the different patterns of 
gene expression. 
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I. Introduction 



Gene expression is the central activity in a living cell. Genes are fragments of DNA 
molecules and provide the blueprint for the synthesis of functional molecules such as 
RNAs and proteins. In each cell, at any instant of time, only a subset of genes present 
is active in directing RNA/protein synthesis. The gene expression is "on" in such a 
case. There are two major steps in gene expression: transcription and translation. 
During transcription the sequence along one of the strands of the DNA molecule is 
copied or transcribed in a RNA molecule (mRNA). During translation, the sequence 
of the mRNA molecule is translated into the sequence of amino acids constituting 
a protein, i.e., a protein molecule is synthesized. Regulation of gene expression is 
an essential process in the living cell and determines the rates and patterns of gene 
expression. An in-depth understanding of gene expression and its regulation is the 
central focus of biology [T]. 

The processes of transcription and translation involve several biochemical reac- 
tions, the kinetics of which determine how the number of participating biomolecules 
changes as a function of time. In the traditional differential rate-equation approach, 
the time evolution of a system of chemical reactions is assumed to be continuous and 
deterministic. In reality, the time evolution is not a continuous process as molecular 
population levels in a reacting system change only by discrete integer amounts. Fur- 
thermore, the time evolution is not deterministic as the collision of molecules which 
brings about a chemical reaction is a probabilistic event. The deterministic rate equa- 
tion approach is justified when the number of molecules of each chemical species is 
large compared to thermal fluctuations in the concentration. In a living cell, the 
number of molecules participating in different biochemical reactions is often small 
and there are considerable fluctuations in the reaction rates. As a result, the time 
evolution of the reacting system, in terms of how the number of reacting molecules 
changes as a function of time is stochastic rather than deterministic. There is now 
an increasing realization that stochasticity plays an important role in determining 
the outcome of biochemical processes in the cell |2J- Stochastic effects in gene ex- 
pression explain the pronounced cell-cell variation observed in isogenic populations. 
A cell may have the option of proceeding along one of two possible developmental 
pathways. The pathway selection is probabilistic and the cell fate depends on the 
particular choice of pathway. Thus, even a clonal population of cells can give rise 
to two distinct subpopulations in the course of time. The randomization of pathway 
choice leads to diversity and increases the likelihood of survival of organisms in widely 
different environments. 

The issue of stochasticity (randomness or noise) and its effect on cellular pro- 
cesses as well as on the operation of synthetic devices like genetic switches, has been 
addressed in several theoretical studies pU IU El El 13 El M EEE HI]- A complete un- 
derstanding of cellular processes requires an appreciation of events at the level of 
an individual cell and subsequent extrapolation to an ensemble of cells. Recent ex- 
perimental advances have made it possible to study processes within a single cell 
unmasked by ensemble averaging [12]. The simplest event one can study at the sin- 
gle cell level is that of the expression of a reporter gene such as lacZ and GFP. In 
the former case, the end product is an enzyme /3-galactosidase which is capable of 
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hydrolyzing a noncolored substrate to a colored product. In the latter case, the pro- 
tein itself is fluorescent. Hence, the gene expression can be directly studied either 
colorimetrically or fluorometrically at the level of an individual cell. Recent experi- 
ments using such techniques, provide evidence that gene expression occurs in abrupt 
stochastic bursts at the level of an individual cell [EH EH EH] • Two very recent experi- 
ments [THE] provide direct evidence of stochasticity in gene expression. In both the 
experiments, a quantitative measure of the noise associated with gene expression has 
been obtained. The noise has both intrinsic and extrinsic components. Intrinsic noise 
is the difference in protein synthesis which arises when two identical copies of a gene 
are expressed under the same conditions. Extrinsic noise occurs due to fluctuations 
in the cellular components required for gene expression. The experiments provide a 
quantitative framework for the characterization of noise in gene regulatory networks. 

Some earlier experiments on both prokaryotic and eukaryotic cells have provided 
evidence of the so-called "all or none" phenomenon in gene expression ITU H5| 
[TBI EH 1213 123] ■ This implies that in an individual cell, gene expression is either 
low/off or has a high value. In an ensemble of cells the protein levels are distributed 
in a bimodal manner, a large fraction of cells synthesize proteins at a low (may be 
zero) level or produce them at a high level. Most of the experiments require the 
presence of inducers/enhancers to observe bimodality. There is strong experimental 
evidence that inducer/enhancers increase the number of expressing cells but not the 
level of expression per cell [21J. The process of gene expression is analogous to a 
binary switch which can be in "on" and "off" positions. Inducer/enhancer molecules 
make it favorable for the switch to be in the "on" position. Some theories have been 
proposed so far to explain the so-called "all or none" phenomenon in prokaryotic gene 
expression. The theories are mostly based on an autocatalytic feedback mechanism 
[IH1 EH EDI 121] , synthesis of the gene product gives rise to the transport or production 
of inducer molecules which in turn promote further gene expression. In section II of 
this paper, we propose a model of gene expression and show that in the presence of 
a sufficient number of inducer molecules in a cell, gene expression in that cell occurs 
at a high level. In cells where inducer molecules are absent or are few in number, 
gene expression occurs at practically zero level. The role of autocatalysis in the 
"all or none" phenomenon is also commented upon. The method employed for the 
study is that of stochastic simulation based on the Gillespie Algorithm (GA) [25] . 
The GA provides a stochastic realization of the temporal pattern of gene expression 
and is more realistic and accurate than the deterministic differential rate equation 
approach. Our model of gene expression includes the major biochemical reactions 
involved in transcription and translation. In section III, we explore the parameter 
space of the model and obtain different temporal patterns of gene expression. One 
parameter region of particular interest corresponds to stochastic flips between the 
states 1 and 2 at random time intervals. In state 1, the protein level is zero. This is 
an example of a binary switch which makes stochastic transitions between the states 
1 and 2 and the temporal process is analogous to a two-state jump phenomenon. The 
effect of changing the reaction parameters on the temporal patterns of gene expression 
is further studied. Some of the results can be understood in the framework of a simple 
mathematical model. 
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II. Stochastic model of gene expression 



We consider a single gene. The gene is transcribed into mRNA by an enzyme called 
RNA polymerase (RNAP) . The process is initiated with the binding of RNAP to a site 
called promoter, usually near the beginning of the transcribed sequence. Expression 
of most genes are regulated at the level of transcription and more specifically during 
the initiation of transcription, that is, before the first phosphodiester bond is formed. 
Regulation of transcription initiation is achieved by the binding of a regulatory protein 
(R) to an overlapping segment of DNA (called operator O) resulting in a turning off of 
mRNA production. RNAP and the regulatory R molecules are mutually exclusive. If 
RNAP binds to the promoter region first, it prevents the binding of R to the operator 
region and vice versa. An inducer molecule (I) may bind to R both when R is free 
and when R is bound to the operator O. In the later case, the complex of I and R 
detaches from the operator. As long as R is forming a bound complex with I, it is 
unable to bind at O and so cannot function as a regulatory protein. R regains its 
activity when the inducer dissociates from the I_R complex and R is able to occupy 
the operator region once more. The biochemical reactions considered in the model of 
gene expression are: 
Reaction 1: 

+ R = 0_R (1) 

Regulatory molecule R binds to the operator region O to form the bound complex 
0_R. 

Reaction 2: 

0_R -> O + R (2) 

Bound complex 0_R dissociates into free R and O. 
Reaction 3: 

P + RNAP = P_RNAP CC (3) 

RNA polymerase (RNAP) binds to the promoter region P forming the closed complex 

P_RNAP CC . 

Reaction 4: 

P_RNAP CC -»• P + RNAP (4) 

The closed complex dissociates into free RNAP and P. 
Reaction 5: 

P_RNAP CC -»• P_RNAP oc (5) 

Isomerization of closed to open complex P_RNAP oc occurs. The open complex is the 
activated form of the RNAP-promoter complex. 
Reaction 6: 

P_RNAP oc -»• TrRNAP + RBS + P (6) 

RNAP clears the promoter region and synthesis of the mRNA chain starts. The 
appearance of the ribosome binding site RBS occurs at the beginning of the mRNA 
chain. TrRNAP denotes transcribing RNA polymerase. 
Reaction 7: 

TrRNAP -> RNAP (7) 
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RNAP completes transcription and is released from DNA. 
Reaction 8: 

RBS + Ribosome -> RibRBS (8) 

Ribosome binds to RBS and RibRBS denotes the bound complex. 
Reaction 9: 

RibRBS -> RBS + Ribosome (9) 

Ribosome dissociates from the bound complex RibRBS. 
Reaction 10: 

RBS — > degradation (10) 

RBS degrades due to the binding of RNAseE at RBS. This binding event is not 
considered separately. 
Reaction 11: 

RibRBS -» EIRib + RBS (11) 

RBS is cleared and the ribosome EIRib initiates translation of mRNA chain. 
Reaction 12: 

EIRib -> protein (12) 

Protein synthesis by transcribing ribosome is completed. 
Reaction 13: 

protein — > degradation (13) 

Degradation of protein product occurs. 
Reaction 14: 

I + R->I_R (14) 

Inducer molecule binds to free regulatory molecule R. I_R is the bound complex of 
I and R. 
Reaction 15: 

LR^I + R (15) 

Bound complex dissociates into free I and R. 
Reaction 16: 

0_R + I^I_R + (16) 

Inducer binds to bound complex 0_R, the complex I_R detaches and the operator 
region O is freed. We emphasize that for many of the steps described above, alternate 
mechanisms exist. However, the mechanism described here is consistent with many 
gene regulatory systems. 

Reaction schemes (1) - (13) are based on those considered in Refs.jU Ej. Tran- 
scription and translation are tightly coupled in prokaryotes. As soon as RNAP leaves 
the promoter region, the 5' end of the mRNA chain, containing the RBS, is available 
for ribosome binding (Reaction 6). This implies that the mRNA chain need not be 
completely synthesized to allow for ribosome binding (Reaction 8), protein synthesis 
by translating ribosome (Reaction 12) and mRNA degradation at RBS (Reaction 10). 
Following Ref.pi], the number of mRNA molecules at any instant of time is given by 
the sum of the numbers of RBS and RibRBS, the bound complex of ribosome and 
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RBS. The functional degradation of mRNA starts at the moment RNAseE binds to 
the RBS. 

We now give a brief description of the Gillespie Algorithm (GA) [2H]. Suppose 
there are N chemical species participating in M chemical reactions. Let X(i), i = 1, 
2, 3 . . . . , N denote the number of molecules of the i th chemical species. Given 
the values of X(i), i = 1, 2, 3, . . . , N at a time t, the GA is designed to answer two 
questions: (1) when will the next reaction occur? and (2) what type of reaction will 
it be? Let the next reaction occur at time t+r. Knowing the type of reaction, one 
can adjust the numbers of participating molecules in accordance with the reaction 
schemes. Thus, with repeated applications of the GA, one can keep track of how the 
numbers X(i)'s change as a function of time due to the occurrence of M different types 
of chemical reactions. Each reaction \i {jj, = 1, 2, . . . , M) has a stochastic rate 
constant c M associated with it. This rate constant has the following interpretation: 

c^dt = probability that a particular combination of reactant molecules participates 
in the [i th reaction in the infinitesimal time interval (t, t + dt). Let /i M be the number 
of distinct molecular combinations for the \x th reaction. Then 

a /1 dt=h ll c fl dt = probability that the \i th reaction occurs in the infinitesimal time 
interval (t, t + dt). 

Let P(r,/i) dr be the probability that the next reaction is of type /i and occurs in the 
time interval (t+ r, t+r + dr). It is straightforward to show that 

P(t, fi) = a M exp(-a r) (17) 

where 

M M 

a = H a " = Yl h " C v ( 18 ) 
u=l v=l 

What is needed now is to generate a random pair (r,/z) according to the probability 
distribution (17). Let ri and r 2 be two random numbers obtained by invoking the 
standard unit interval uniform random number generator. One can then show that r 
and /i are obtained as 

r=(i)ln(i) (19) 
and /i is taken to be the integer for which 

A 4 

a v < r 2fl0 < &v (20) 

A rigorous proof of the formulae (19) - (20) is given in Ref.[25j. Once r and fx 
are known, the time evolution of the reacting system is specified. The stochastic rate 
constant c M is related to the more familiar deterministic reaction rate constant 
through simple relations. In the case of first order reactions, both constants have the 
same value. In the case of second order reactions, the rate constant is divided by 
the volume of the system. We have applied the GA to our model of gene expression 
involving N = 15 types of biomolecules participating in M = 16 biochemical reactions. 
In the initial state (time t=0), the number of free operator and promoter sites is 1. 
The number of R, I, RNAP, and ribosome molecules is Nr = 20, iVj = 20, Nrnap = 
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400 and Nrh, = 350 respectively. The number of all the other biomolecules is set to 
zero at t = 0. The simulation time is up to 2000s, i.e., less than the cell generation 
time typically in the range 2000-3000s. As already mentioned, knowledge of r and 
// enables one to calculate the numbers of biomolecules at time t+r and in this way, 
through repeated applications of the GA, one can keep track of how the different 
numbers change as a function of time. Fig.l shows the number of protein molecules 
present in the system as a function of time. The stochastic rate constants of the 
sixteen reactions are Ci = 0.5, C2 = 0.004, C3 = 0.02, C4 = 0.001, C5 = 0.8, C6 = 0.9, 
c 7 = 0.08, c 8 = 0.01, c 9 = 0.001, cio = 0.3, c u = 0.8, c 12 = 0.7, C13 = 0.003, c 14 = 
0.7, C15 = 0.001 and Ci6 = 0.3 respectively. The topmost curve corresponding to Ni 
shows that the protein number reaches a steady level with fluctuations around the 
mean. The dashed curve at the bottom of the figure shows the number of proteins as 
a function of time when the number of inducer molecules Nj (=3) is much less than 
Nr. The protein number in this case is very small, practically zero. The other curves 
corresponds to values of Nj closer to Nr. As Ni approaches Nr, protein levels are 
intermediate between high and low. Figure 2 shows the number of protein molecules 
as a function of time in a different parameter region. 

1800 1 1 1 1 1 1 1 1 1 1 1 




Time in sec 



FIG. 1. No. of protein molecules as a function of time. The stochastic rate constants 
are c x = 0.5, c 2 = 0.004, c 3 = 0.02, c 4 = 0.001, c 5 = 0.8, c 6 = 0.9, c 7 = 0.08, 
c 8 = 0.01, c 9 = 0.001, do = 0.3, di = 0.8, c 12 = 0.7, c 13 = 0.003, c u = 0.7, 
d 5 = 0.001, de = 0.3; N R = 20, N RNAP = 400 and N Rib = 350. 

The "all or none" phenomenon has been observed in both prokaryotic and eukary- 
otic cells. In prokaryotes, genes are often arranged in operons,i.e, sets of contiguous 
genes which include structural and regulatory sequences. A well-known example is 
that of the E. coli lactose (lac) operon [T]. Lac operon consists of three structural 
genes (z, y and a) which code for the three proteins: /3-galactosidase, the enzyme 
that catalyzes the hydrolysis of lactose to glucose and galactose; permease, a carrier 
protein responsible for membrane transport of lactose into the cell and a third protein 
transacetylase. The lac operon contains three regulatory sequences, i, P and O, which 
control the transcription of mRNA leading to the synthesis of the three proteins. The 
sequence i corresponds to a gene lad which is transcribed continuously to synthesize 
a repressor protein, lac repressor, at a low level. Lac repressor binds to the operator 
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FIG. 2. No. of protein molecules as a function of time. The stochastic rate constants 
are the same as in Fig. 1 except that c\ = 0.35 and cu = 0.3. 

sequence O and prevents the transcription of the genes z, y and a so that the /3- 
galactosidase enzyme and the permease molecules are not produced. If the bacterium 
is to grow on lactose (milk sugar) which acts as its carbon source, /3-galactosidase 
must be made available to split the sugar into glucose and galactose. The breakdown 
product of lactose act as an inducer molecule. The inducer attaches to the repressor 
molecule, causing it to release the DNA so that transcription of the structural gene is 
possible. The repressor is freed of the inducer when the lactose supply is exhausted 
and switches off the expression of the structural genes once more. Our simple model 
of gene expression incorporates some of the key features of the lac operon. The role 
of the lac repressor is played by the regulatory molecule R though its synthesis is not 
explicitly considered. There is a single gene in our model analogous to the structural 
gene z expressing the enzyme /3-galactosidase. The inducer molecule I acts in the 
same manner as in the case of the lac operon. 

The experimental observation of the "all-or-none" phenomenon in the lac operon 
has been attributed to autocatalytic feedback mechanisms [TBI EH EI| • At low inducer 
concentrations, some of the bacterial cells synthesize protein at the full rate whereas 
the other cells are in the "off" state. When inducer is added to the colony of bacterial 
cells, simultaneous production of the /3-galactosidase enzyme and permease molecules 
occurs. The permease molecules transport lactose into the cell raising the internal 
inducer concentration which in turn promotes the production of more /3-galactosidase 
and permease. Thus an autocatalytic feedback process is at work and within a short 
time after the appearance of the first permease molecules, the bacterial cell becomes 
fully induced synthesizing the structural proteins at maximum rate. Siegele and Hu 
[20J carried out experiments on gene expression from plasmids containing the araBAD 
promoter in the presence of subsaturating concentrations of the inducer arabinose. 
Again, as in the case of the lac operon, it has been suggested that an autocatalytic 
induction mechanism, due to the accumulation of inducer molecules by transport, is 
at work. However in all the experiments, full induction of cells has been observed even 
in the absence of autocatalysis, i.e., when the inducer availability is not linked to that 
of synthesized protein molecules like permease. Our model of gene expression does 
not include an autocatalytic feedback process and a detailed analysis of the simulation 
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results shows that if the number of inducer (I) molecules, iVj, is greater than or equal 
to the number of regulatory (R) molecules, Nr, in a cell, the cell reaches a steady 
state which is state 2 (high protein level). If Nj= or <C Nr, the cell is in state 1 
(low/zero level) but as Nj approaches Nr , protein levels intermediate between high 
and low are obtained. The "all or none" phenomenon becomes more pronounced in 
the presence of the autocatalytic induction (AI) mechanism. If Nj is originally small 
in a cell, the AI mechanism leads to a rapid increase in Nj . When Nj is > Nr, the 
cell is in state 2 in the steady state. The magnitude of the protein level in state 2 is 
independent of the value of Nj in agreement with experimental results. Autocatalysis 
is responsible for a sharp bimodal distribution in protein levels. In the absence of 
autocatalysis, the magnitude of the protein level in state 2 is the same as in the case 
of autocatalysis but the bimodal distribution become less sharp due to intermediate 
protein levels in a fraction of cells. Our simulation results are in agreement with 
experimental observations |THl EH EH ■ Experiments [131 EH EH] on single mammalian 
cells (eukaryotic cells) have provided further evidence of bimodality in the distribution 
of protein levels in an ensemble of cells. Again, one finds that the amount of enhancer 
affects the number of expressing cells but not the level of expression. In other words, 
the enhancer increases the probability rather than the rate of transcription. In the 
experiments carried out by Zlokarnik et al. the reporter gene synthesizes the 

protein /3-galactosidase. In unstimulared cells, the number of these proteins is low, 
in the range 150 - 300. Under the action of the stimulating agent carbachol, rapid 
conversion to a state, in which 15000 - 20000 /3-lactonase molecules are present in 
a cell, is obtained. The major conclusion of the single cell experiments mentioned 
above is that in the systems considered, the cellular state is bistable. A cell can 
exist in two stable steady states: gene expression "off/low" and gene expression "on" 
with a high level of protein production. Addition of inducer/enhancer to the system 
increases the fraction of cells in the "high" state. In the case of eukaryotic systems, 
however, the mechanism of enhancer action is not well understood. Enhancers have 
been suggested to give rise to two major types of response. Enhancers increase the 
rate of transcription through mainly enhancing the rate of close to open complex 
formation of RNAP bound to the promoter region [2E1- The second type of response 
is of the "all or none" type (211 EH]. Enhancers in this case do not increase the rate 
of transcription but increase the fraction of cells in the high state (state 2). The "all 
or none" phenomenon observed in some eukaryotic systems [T2l EH E3 [21] dose not 
involve autocatalysis explicitly. Thus, a more general mechanism than in the case of 
prokaryotic systems is required to explain the bimodal distributions in protein levels 
in an ensemble of cells. 



III. Patterns of gene expression 

We now consider the model of gene expression in the absence of inducer molecules. 
The values of M, the total number of reactions and N, the number of different types 
of biomolecules, are both thirteen. To explore the full parameter space, one has 
to treat the thirteen stochastic rate constants corresponding to the same number 
of reactions as variables. Experimental results wherever they are available, show 
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that the usual rate constants, k^'s, to which the stochastic rate constants c M 's are 
related, can vary over a wide range depending on the type of gene and the nature of 
the cellular environment [2H]. Since the exploration of the full parameter space is a 
daunting task, we report on some of the more general patterns of gene expression in 
a restricted subspace. The effect of changing the stochastic rate constants on specific 
patterns is also studied. We have included an optional feature in our model of gene 
expression, namely, cooperative binding of RNAP to the promoter region P. The 
possibility of such a binding has been suggested earlier in a simplified probabilistic 
model of gene expression [2Z|- In the present model, cooperative binding implies that 
the rate constant for the binding of RNAP at P is enhanced by a factor q if the 
binding event is immediately preceded by the Reaction 6 in which a RNAP clears the 
promoter region. Cooperative binding of proteins to DNA is now well established. 
In most cases of regulatory proteins, the binding cooperativity is mediated through 
protein-protein interaction although increasing evidence of DNA mediated effects have 
been reported [2H|- Some recent experiments on the prokaryotic system E. coli have 
shown that the transcriptional activity of the promoter is intrinsically sensitive to 
the superhelical density of the DNA template [21112011211122]. In fact as pointed out 
by McClure[26j, supercoiling gives rise to considerably more diversity in the patterns 
of promoter strength (the ability to bind weakly or strongly) than do mutations of 
auxiliary proteins. There is now experimental evidence that transcription generates 
increased negative supercoiling through several hundred base pairs[29, 30, EH [321 
133] . In principle, this can facilitate the binding of RNAP to the DNA or decrease 
the energy of activation required for the isomerization of RNAP-promoter complex 
from closed to open formp3|. Thus, it is entirely plausible and likely that active 
transcription downstream of the promoter site may lead to increased binding of RNAP 
(cooperative binding) and enhanced rate of open complex formation (stochastic rate 
constants c 5 ). 

We now describe three major types of gene expression pattern as a function of 
time. As before, the GA is applied to the system of thirteen biochemical reactions 
constituting the processes of transcription, translation and the regulation thereof. In 
the gene expression pattern designated as Type A, the protein production occurs in 
abrupt stochastic bursts. A variable number of proteins is produced in each burst. 
The type A pattern of gene expression has been observed experimentally^ [TJJ [TH] 
and has been attributed to stochastic effects. Figures 3 and 4 show two such patterns 
of protein production for q = 1 (no cooperative binding). 

In the Type B pattern, the protein level reaches a steady state with fluctuations 
around the mean (Figures 5 and 6). This type pattern is quite common and is 
routinely observed in experiments. 

The Type C pattern of gene expression has an interesting structure. As in the 
case of the Type A pattern, protein production occurs in stochastic bursts, i.e., at 
random time intervals. The bursts may be of various durations but in each burst, the 
protein number attains the same level (with attendant fluctuations) in a very short 
time. Similarly, the decay of the protein level from high to zero occurs in a small time 
interval. Figures 7(a) and 7(b) show the patterns of mRNA and protein production 
as a function of time. The number of regulatory (R) molecules (Nr), RNAP (Nrmap) 
and ribosome (Nrm) is 10, 400 and 200 respectively. 
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FIG. 3. No. of protein molecules as a function of time. The stochastic rate constants 
are a = 0.008, c 2 = 0.004, c 3 = 0.007, c 4 = 0.001, c 5 = c 6 = 1, c 7 = 0.4, c 8 = 0.01, 
c 9 = 0.001, cio = 0.1, c n = ci 2 = 1, ci 3 = 0.03; iV R7VA p = 400, = 10 and 
^ = 200. 




FIG. 4. No of protein molecules as a function of time. The stochastic rate constants 
are a = 0.008, c 2 = 0.004, c 3 = 0.08, c 4 = 0.001, c 5 = c 6 = 1, c 7 = 0.4, c 8 = 0.01, 
c 9 = 0.001, cio = 0.3, c u = c 12 = 1, c 13 = 0.05; Nrnap = 400, JV fl = 10 and 
N Rib = 200. 
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FIG. 5. No. of protein molecules as a function of time. The stochastic rate constants 
are a = 0.008, c 2 = 0.004, c 3 = 0.5, c 4 = 0.001, c 5 = c 6 = 1, c 7 = 0.4, c 8 = 0.01, 
c 9 = 0.001, cio = 0.01, c u = cia = 1, ci 3 = 0.005; N RNAP = 400, N R = 10 and 
iV^ = 200. 
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FIG. 6. No. of protein molecules as a function of time. The stochastic rate constants 
and the other parameter values are the same as in Fig. 5 except that cio = 0.1, and 
C13 = 0.002. 
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FIG. 7. No. of mRNA (a) and protein (b) molecules as a function of time. The 
stochastic rate constants are C\ = 0.01, c 2 = 0.004, c 3 = 0.7, c 4 = 0.001, c 5 = c 6 = 20, 
c 7 = 0.4, c 8 = 0.01, c 9 = 0.001, cio = 0.2, c n = c 12 = 1, c 13 = 0.08; N RNAP = 400, 
N R = 10 and N Rlb = 200. 
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We now discuss the physical origin of the Type C pattern. The stochastic rate 
constant C3 for RNAP binding (Reaction 3) is considerably higher than that of the 
binding of R at the operator site O (Reaction 1). The initial number of RNAP 
molecules is also larger than that of R molecules. As explained in section II, a^dt 
= h^c^dt is the probability that the \x th reaction occurs in the infinitesimal time 
interval (t , t + dt ) . The number of distinct molecular combinations for the fi 
th reaction is 10 and 400 for the Reactions 1 and 3 respectively. The corresponding 
stochastic rate constants have the values ci= 0.01 and C3 = 0.7. Thus, Reaction 3 is 
more probable than Reaction 1. Note that the stochastic rate constants C5 , C6 are 
considerably high. Reactions 5 and 6 are associated with the transcription process, 
namely, isomerization of the closed complex of RNAP bound to the promoter region 
P to the open complex and subsequent clearance of the promoter region by RNAP. 
After the binding of a RNAP to P, a host of factors including high values of some 
of the relevant rate constants, leads to a sharp rise in the number of proteins to a 
level determined by the transcription, translation and protein degradation rates. The 
protein level is maintained over a time interval due to the balancing of the rates of 
synthesis and degradation. Binding of the R molecule to O, though less probable 
than that of RNAP at P, can occur with a finite probability. Once the R molecule is 
bound to O, it continues to remain bound for some time as the dissociation rate (02 
= 0.004) is low. This prevents the binding of a RNAP to P during the time interval 
in which R stays bound to O, leading to a sharp fall in the number of proteins to 
zero. The subsequent dissociation of the R molecule from the operator O, followed by 
the binding of a RNAP to the promoter region P, tilts the balance in favour of state 
2. The cellular state thus flips from state 1 (no. of proteins zero) to state 2 (no. of 
proteins high) and vice versa at random time intervals dictated by stochastic binding 
and dissociation events at O. In section II, we have discussed a bimodal distribution 
in protein levels due to the "all or none" phenomenon in an ensemble of cells. In 
the present bimodal distribution in protein levels is obtained in an ensemble 

of cells if protein levels are measured at particular instant of time. The heights and 
widths of the two peaks may change as a function of time but the distribution remains 
bimodal. 

400 
350 
300 
B 250 
200 
150 
100 
50 


0.2 0.4 0.6 0.8 1 1.2 

m 

FIG. 8. Distribution of the no. N(m) of cells expressing fraction m of the average 
number of proteins. The total number of cells is 1000. The stochastic rate constants 
are the same as in Fig. 7. 
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FIG. 9. No. of protein molecules as a function of time. The stochastic rate constants 
are a = 0.008, c 2 = 0.004, c 3 = 0.04, c 4 = 0.001, c 5 = c 6 = 1, c 7 = 0.4, c 8 = 0.01, 
c 9 = 0.001, cio = 0.1, c u = c 12 = 1, c 13 = 0.03; N RNAP = 400, N R = 10 and 
N Rib = 200. 

Figure 8 shows the distribution of the number N(m) of cells expressing a fraction 
m of the average number of proteins at a particular instant of time. The total number 
of cells is 1000. The values of the stochastic rate constants are the same as in Figs. 
7(a) and 7(b). Type C pattern of gene expression is also obtained for lower values of 
c 5 and c 6 (Fig. 9) though better quality patterns are obtained for high values of the 
rate constants. Figures 10(a) and 10(b) show the temporal variations of the mRNA 
and protein numbers with stochastic rate constants the same as in Figs. 7(a) and 
7(b) but the enhancement factor q, associated with cooperative RNAP binding, has 
been raised from 1 to 10. Comparing the two sets of Figures, one concludes that 
cooperativity increases the duration of state 2 ("high" level). The magnitude of the 
mean level remains unchanged. In an ensemble of cells, a greater fraction of cells is 
in state 2 than in the earlier case. The origin and nature of bimodal distribution in 
protein levels are different for the model system considered in section II and the type 
C pattern of gene expression. In the latter case, no autocatalytic feedback is necessary 
to obtain a bimodal distribution. Positive (autocatalytic) feedback mechanism has 
been invoked to explain the "all-or-none" phenomenon in prokarvoticfTHt fT9l 20j and 
eukaryotic[2ni systems. Experimental reports of the phenomenon in some eukaryotic 
systems [THl CU EH] suggest that autocatalytic (positive) feedback is not essential for 
the obsevance of the phenomenon. In these systems, activation to the high protein 
level is enhancer mediated. The origin of bimodal distribution in protein levels is, 
however, yet to be elucidated. Figure 11 shows the number of proteins (solid line) 
and mRNA molecules (dotted line) as a function of time to make a simultaneous 
comparison of the rise and decay of protein and mRNA levels. Note that after the 
number of mRNA molecules become zero, there is a time delay before the number of 
proteins falls to zero. Even when there are no mRNA molecules in the system, some 
proteins remain which degrade to zero level at a decay rate lower than that of the 
mRNA molecules. In Fig. 12, a pattern of gene expression is shown in which the 
protein number never falls to zero and a variable number of proteins is synthesized 
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FIG. 10. No. of mRNA (a) and protein (b) molecules as a function of time. The 
stochastic rate constants are c\ = 0.01, c 2 = 0.004, c 3 = 0.7, c 4 = 0.001, c 5 = c 6 = 20, 
c 7 = 0.4, c 8 = 0.01, c 9 = 0.001, cio = 0.2, c n = c 12 = 1, c 13 = 0.08; N RNAP = 400, 
N R = 10, N Rib = 200 and q = 10. 

as a function of time. 

We have further checked how robust pattern C ( the gene expression pattern shown 
in Figure 7) is when the different stochastic parameters are changed. 




FIG. 11. No. of protein and mRNA molecules as a function of time. The stochastic 
rate constants are C\ = 0.01, c 2 = 0.004, c 3 = 0.5, c 4 = 0.001, c 5 = c 6 = 20, c 7 = 0.4, 
c 8 = 0.01, c 9 = 0.001, cio = 0.2, c n = c 12 = 1, c 13 = 0.08; A^ap = 400, N R = 10, 
A ffife = 200. 

We change one parameter at a time keeping all the other parameter values the 
same as in Fig. 7. As ci decreases from 0.01 (ci = 0.01 in Fig. 7), the duration of 
state 2 increases and ultimately state 2 becomes the steady state (Figs. 13(a) and 
13(b)), i.e., a Type B pattern is obtained. As Ci increases from 0.01, the Type C 
pattern is still obtained (Fig. 13(c), ci = 0.02) but for higher ci values, say, ci = 
0.03, the pattern of gene expression becomes of Type A (Fig. 13(d)). If c 2 is varied, 
then as c 2 decreases from 0.004 (Fig. 7), the total duration of state 2, for t in the 
range - 2000s, decreases (Fig. 13(e) with c 2 = 0.008). If c 4 is varied, one finds that 
the Type C pattern of gene expression is obtained over a wide range of values. The 
greater the values of C5 and C6 , the higher is the level of proteins attained in state 
2. The best Type C patterns are obtained if the values of C5 and C6 are high. Figs. 
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FIG. 12. No. of protein molecules as a function of time. The stochastic rate constants 
are a = 0.008, c 2 = 0.004, c 3 = 0.08, c 4 = 0.001, c 5 = c 6 = 1, c 7 = 0.4, c 8 = 0.01, 
c 9 = 0.001, cio = 0.1, cn = cia = 1, cia = 0.03; N RNAP = 400, Afo = 10, N mb = 200. 

13(f) and 13(g) correspond to c 5 = 2, c 6 = 20 and c 5 = 20, c 6 = 2 respectively. If C7 
is varied, the Type C pattern is obtained over a wide range of values. If c 8 increases 
from the value 0.01 (Fig. 7), Type C pattern becomes Type A (Fig. 13(h), c 8 = 
0.015). If c 8 decreases from the value 0.01, the total duration of state 2 increases(Fig. 
13(i)). If cio decreases from the value 0.2 (Fig. 7), the Type C pattern is lost (Fig. 
1 3 ( j ) , Cio = 0.1). If Cio increases from 0.2 the Type C pattern is obtained over a wide 
range of values of Ci . The best type of Type C pattern is obtained for similar values 
of Cn and C12. 
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FIG. 13(a). No. of protein molecules as 
a function of time. The stochastic rate 
constants are C\ = 0.002, c 2 = 0.004, c 3 = 
0.7, c 4 = 0.001, c 5 = c 6 = 20, c 7 = 0.4, 
c 8 = 0.01, c 9 = 0.001, cio = 0.2, cn = 
C12 = 1, cia = 0.08; Nrnap = 400, AT* = 
10 and iV ffife = 200. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that c\ = 
0.002. 
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FIG. 13(b). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 13(a) except that 
ci = 0.0008. 
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FIG. 13(c). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 13(a) except that 
ci = 0.02. 
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FIG. 13(d). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 13(a) except that 
ci = 0.03. 
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FIG. 13(e). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that c 2 = 
0.008. 
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FIG. 13(f). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that C5 = 
2 and c 6 = 20. 
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FIG. 13(g). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 13(f) except that 
c 5 = 20 and c 6 = 2. 
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FIG. 13(h). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that c§ = 
0.015. 




200 400 600 800 1000 1200 1400 1600 1800 2000 



Time in sec 

FIG. 13 (i) . No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that c 8 = 
0.004. 
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FIG. 13(j). No. of protein molecules as 
a function of time. The stochastic rate 
constants and the other parameter values 
are the same as in Fig. 7 except that c w = 
0.1. 
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We have also studied the effect of changing the number Nr, Nr^ap and Nro, of 
R, RNAP and ribosome molecules respectively. As Nr increases from 10 (Fig. 7), the 
Type C pattern is gradually lost. Reducing Nr increases the total duration of state 
2. lfN R = 0, state 2 becomes the steady state. If N RN ap is reduced from 400 (Fig. 7) 
the Type C pattern worsens gradually. If N RN ap increases beyond 400, the Type C 
pattern is still obtained. The magnitude of the protein level in state 2 remains more 
or less the same. If N Rib decreases from 200, the total duration of state 2 increases. 

Some general conclusions that can be made, on the basis of the stochastic simula- 
tion of gene expression patterns, are as follows. The protein level attained in state 2 
is not affected by changes in the rate constants ci, C2, C4, C7, en, C12 and the enhance- 
ment factor q. The rate constants C5, eg, eg, C10 and C13 determine the magnitude of 
the protein level. If ci decreases and C2 increases, the total duration T of state 2 in 
the time interval 0-2000s increases. Transition from Type C — >Type B can occur by 
decreasing ci(Fig. 13(b)), increasing q and by making the values of C5 and C6 unequal 
(Fig. 13(g)). Transition from Type C — ► Type A pattern can occur by increasing 
Ci (Fig. 13(d)), c 8 (Fig. 13(h)) or by decreasing Ci (Fig.l3(j)). Type C patterns are 
favourable for high values of C5 and C6 with cs~ C6 and also en — C12. 

Kepler and ElstonfHj have proposed a simple mathematical model of gene expres- 
sion with no feedback. In their model, if an activator molecule occupies the operator 
region, protein production occurs at a rate ot\. If the operator is unoccupied by the 
activator molecule, protein production occurs at a lower rate ckq. The model does 
not include the intermediate steps of protein synthesis like transcription, ribosome 
binding for the initiation of translation etc. The model, however, provides lots of 
physical insight on stochastic effects in the form of fluctuations in the discrete states 
of the operator (unoccupied/occupied) on gene expression. Chemical reactions that 
change the state of the operator are termed operator fluctuations. Approximations to 
the dynamics were made for the cases in which the protein number is large or the 
operator fluctuations are fast. In the first case, the effective rate of protein synthesis 
fluctuates randomly in time between high (synthesis rate ot\) and low (synthesis rate 
«o) levels. In the latter case, the fluctuations are effectively averaged out over larger 
time scales. In our model, the major biochemical reactions/events in protein synthe- 
sis have been included and instead of activators we have regulatory molecules which 
act as repressors. The simulation based on GA provides an accurate knowledge of 
the microscopic origins of the different types of temporal gene expression patterns. 
A simple mathematical model provides insight on the origin of different patterns and 
transition from one type pattern to another. Let m denote the number of proteins at 
time t devided by the maximum number of proteins. The equation 

dm 
~~dt 

describes the rate of change in the number of proteins. The possible values of x are 
and 1 so that in the steady state m can be either 1 (high level, state 2) or (low 
level, state 1). The variable x randomly switches between the two states and the 
transition rate from state 1 — > state 2 is T\ and that from state 2 — > state 1 is r 2 . 
Let Pj(m = M,t), j = 1,2, be the probabilities of being in the state j. One can then 
write down the master equations[34j 
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(21) 



(22) 



dPi d((l - M)P 1 ) 



- r 2 P 1 + nP 



dt dt 

The steady state distribution of P(m,t) = P (m,t) + Pi(m,t) is given by 



P(m — M) — AM Tl ~ (I - M) 



\T2-1 



(23) 



(24) 



where A is the normalization constant. Figs. 14(a), (b) and (c) show the distribution 
P(m) versus m for different values of r\ and r 2 . The distribution in Fig. 14(a) 
corresponds to the Type C pattern of gene expression (the magnitude of the high 
protein level in state 2 is normalized to the value 1) in the steady state and is bimodal 
in nature. Fig. 14(b) describes the Type B pattern of gene expression and there 
is a single peak corresponding to the steady state level, the magnitude of which is 
normalized to 1. Fig. 14(c) shows a broad distribution in protein levels. Fig. 12 shows 
a pattern of gene expression based on simulation of the detailed model which gives 
rise to a probability distribution of protein levels similar to that shown in Fig. 14(c). 
The transition rates r\ and r 2 are functions of the different stochastic rate constants 
though the actual functional relationship is yet to be worked out. The rate r\ depends 
dominantly on c 2 , C5, Cq, c§ and C10 whereas the rate r 2 is mainly determined by C\ 
and C13. The simple model illustrates how the different probability distributions arise 
and the transition from one type to another occurs as the transition rates r\ and 
r 2 are changed. For slow transition rates, the states 1 and 2 can be distinguished 
and the probability distribution P{m) is bimodal. For fast transition rates, m has 
values intermediate between and 1 and the bimodality is smeared into a single 
broad distribution. For large r\ and small r 2 , the peak is around the high value 1. 
The transitions from one type of probability distribution to another as r\ and r 2 are 
varied are consistent with the changes in the nature of the temporal patterns of gene 
expression (Figs. 12 and 13), brought about by changes in the various stochastic rate 
constants on which the transition rates r\ and r 2 of the simpler model depend. 
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FIG. 14(a). Distribution of P{m) as a function of m for r\ = 0.1, and r 2 = 0.01. 
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FIG. 14(b). Distribution of P{m) as a function of m for r\ = 2.0, and r 2 = 0.2. 
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FIG. 14(c). Distribution of P{m) as a function of m for r 1 = 1.3, and r 2 = 1.01. 

IV. Concluding remarks 

In this paper, we have considered two models of gene expression in a single cell. 
In the first model, the effect of inducer molecules is explicitly considered and the 
cellular kinetics are described by a set of sixteen biochemical reactions. In the second 
model, the cellular kinetics are described by a set of thirteen biochemical reactions. 
Both the models include the basic steps of transcription and translation as well as 
transcriptional regulation through the binding of a regulatory molecule R to the 
operator region, with R acting as a repressor of transcription. Additionally, in the 
second model a transcriptionally generated cooperativity factor q in the binding of 
RNAP has been introduced as an optional feature. The patterns of gene expression as 
a function of time have been determined in both the models with the highly accurate 
stochastic simulation method based on the Gillespie Algorithm. Analytic approaches 
[HI E3 GDI EH are possible in studying the temporal evolution but in these cases the 
detailed biochemical reactions are lumped together into a few effective processes. This 
makes it possible to determine the temporal evolution of the system in a chemical 
Master Equation approach which is stochastic in nature. Alternatively, the effect 
of stochasticity can be taken into account by the inclusion of "noise" terms in the 
differential rate equations. The GA provides a more detailed picture of the kinetics 
though the computational efforts required in its implementation are considerable. In 
section II, we have studied gene expression in the presence of inducer molecules. The 
major result of this study is that when the number of inducer molecules, iVj, is greater 
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than or equal to the number of regulatory molecules, Nr, the cellular steady state 
is state 2 in which the proteins are synthesized at a high level. When Ni <C Nr , 
the cellular state is state 1 in which the protein level is low/zero. As Nj approaches 
Nr, protein production occurs at intermediate levels. Fig.l provides an example of 
how drastic the effect of even a single molecule can be on gene expression. There 
is a considerable difference of protein levels attained in the steady state when Nj is 
changed from 19 to 20. This single molecule effect can be observed over a wide region 
of parameter space and may give rise to threshold phenomena. These results are new 
and have not been reported earlier. McAdams and Arkin[2j have pointed out that 
even a single molecule can switch the biochemical state of a cell. Togashi and Kaneko 
[35j have studied an autocatalytic reaction system with a small number of molecules 
and shown that due to the nonlinear dynamics, amplification of small changes can 
give rise to single molecule switches. 

The results of our simulation explain why an autocatalytic induction mechanism 
gives rise to the "all-or-none" phenomenon observed in cells in the presence of inducer 
molecules[l8l fT9| l2fl]. Due to autocatalytic induction, the small number of inducer 
molecules initially present in a cell is quickly amplified so that the condition Nj > Nr 
is satisfied in a short interval of time. The cell then exists in state 2 with a high protein 
level. Other cells in which iVj is zero exist in state 1 with zero protein level. In the 
presence of a subsaturating concentration of inducer molecules, the distribution of 
protein levels in an ensemble of cells is bimodal. In the absence of autocatalytic 
induction, protein production occurs at low, high as well as at intermediate levels 
so that the bimodal distribution gets smeared. This is in keeping with experimental 
observations. 

The model studied in section III does not include inducer molecules. Simulation 
based on the GA shows the existence of three types of pattern of gene expression 
as a function of time. In the Type A pattern, protein synthesis occurs in abrupt 
stochastic bursts and a variable number of protein is produced in each burst. There 
is considerable experimental evidence jTHl EH EH] for this type of gene expression and 
in section III, some examples of this type of pattern have been given (Figs. 3 and 4). 
In Type B pattern, protein levels reach a steady state (Figs. 5 and 6). This is very 
common type of pattern observed in gene expression experiments. 

In the Type C pattern (Fig. 7), the cellular state makes random transitions be- 
tween states 1 and 2 . This is a different manifestation of the "all or none" phenomenon 
and measurement levels in an ensemble of cells is bimodal (Fig. 8). Type C patterns 
of gene expression have been obtained by Kepler and Elston |H] in their study of 
model gene expression systems (see Fig. 4 of Ref. |Hj) using the Master Equation 
Approach. In these models detailed biochemical reactions have been replaced by a 
few effective processes so that the models are mathematically tractable. Kepler et 
al. have specifically considered the effect of fluctuations in the state of the opera- 
tor on gene expression. Their conclusion is that the operator fluctuations can induce 
bistability in parameter regions which gives rise to monostability in the deterministic, 
i.e., the zero noise limit or destroy bistability if it exists in the noise- free case. In 
the deterministically bistable region, the gene acts like a genetic switch and external 
noise/perturbation is needed to flip the switch from one state to the other. This is the 
principle of operation behind the noise-based switches and amplifiers for gene expres- 
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sion proposed by Hasty et al [TT]. If a system is stochstically bistable, the fluctuations 
in the system flip the switch between the two states (say, state 1 and state 2) at ran- 
dom time intervals. As in Ref. (Hj, the operator fluctuations have been explicitly 
considered in the models studied in Sections II and III and it has been shown that in 
certain parameter regimes, stochastically bistable behaviour corresponding to Type 
C pattern of gene expression is obtained. Our stochastic simulation results, obtained 
by using the highly accurate Gillespie Algorithm, provide a microscopic basis as well 
as quantitative estimates of the different rate constants for obtaining Type C pattern 
of gene expression. The simulation method can be used to study a large number of 
reactions which is not feasible in the formalism of a mathematical model. The signif- 
icant omission in Ref. [Hj is that no distinction has been made between transcription 
and translation. The models describe direct translation from gene into protein. As 
pointed out in the paper, the simplification may have considerable impact on cellular 
phenomena. For example, the synthetic repressilator network [6J would not oscillate 
if a time delay between transcription and translation were absent. In our model of 
gene expression, all the major biochemical reactions involved in transcription and 
translation have been explicitly taken into account and the time delay is clearly seen 
in Fig. 11. There is some experimental evidence [23 EDI EU E2J of transcriptionally 
generated cooperative binding of RNAP to the promoter region. For non-zero valus of 
the cooperativity factor q in our model the duration of state two (high protein level) 
is found to increase. The result is new and the effect of cooperative RNAP binding 
on gene expression needs to be investigated in greater detail. 

We have further studied the effect of changing the various stochastic rate con- 
stants on the temporal pattern of gene expression. Transitions involving Type C — ► 
Type B and Type C — > Type A patterns of gene expression have been obtained by 
changing appropriate stochastic rate constants. These transitions can be understood 
in the framework of a simple mathematical model which does not include the detailed 
biochemical reactions. The model describes a simple gene expression system in which 
random transitions occur between states 1 and 2, corresponding to low and high lev- 
els of protein production. The model has two parameters t\ and r 2 which are the 
transition rates from state 1 — > state 2 and state 2 — > state 1 respectively. In the 
steady state, the probability distribution of protein levels can be calculated. A bi- 
modal distribution corresponds to Type C patterns whereas an unimodal distribution 
is obtained for Type B patterns. A broad distribution in protein levels (Fig. 14(c)) 
is obtained corresponding to the gene expression pattern shown in Fig. 12. Transi- 
tion from one type of distribution to another can be obtained by changing the rate 
constants T\ and r 2 . While experimental evidence for Type A and Type B patterns 
of gene expression is considerable, we do not know of specific experiments exhibiting 
Type C patterns. As mentioned in the Introduction, the Type C pattern is similar 
to that obtained in the case of two-state jump processes in which transitions between 
two states occur at random time intervals (33J EE]- An example is provided by a 
spin- 1 in the presence of a magnetic field and in contact with a heat bath (HE!- The 
literature on two state jump processes is large and one measurable quantity of inter- 
est is the mean first-passage time (MFPT). The time required to switch between two 
states is a random variable and is known as the first first passage time. Determina- 
tion of the MFPT and other characteristic measures of the two-state jump processes 
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describing Type C gene expression has not been attempted in this paper. The Type 
C pattern is reminiscent of a binary digital pulse with states 1 and 2 corresponding 
to the "0" (OFF) and "1" (ON) states. Genes with expression pattern of Type C 
may be combined together to construct binary logical circuits. Recently, synthetic 
networks of genes displaying features of binary logic circuits have been constructed 
[37] . Bialek [113] has studied stability and noise in biochemical switches and has shown 
that switches with long periods of stability and switchability in milliseconds can be 
constructed from fewer than a hundred molecules. The conclusion is arrived at by 
studying a model of the synthesis of a single biochemical species, say, proteins in the 
Langevin formalism. The result obtained is of considerable interest but needs to be 
verified in a detailed approach involving intermediate processes. The correlation of 
the amount of random variation in protein distribution, measured by the Fano factor, 
with the transcriptional and translational rates can be determined in the stochastic 
simulation approach and the results compared with those obtained in experiments 
[THj IT7] . A detailed study on the dominant contributions to noise in our models of 
gene expression is in progress and the results will be reported elsewhere. 
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